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I. INTRODUCTION 



Breakthroughs in numerical relativity during this decade have made it possible to sim- 
ulate, via evolution of the full 3D Einstein equations, binary black hole dynamics through 
inspiral, merger and ringdown of the remnant sing le black hole [MEI (see e. g. recent 



reviews [16l . Il7j). Inspiraling binaries are among the most promising sources of g ravi- 



tational waves for the network of laser interferometric detectors such as LIGO 
VIRGO 



18 and 



19, 20 



Through the construction of templates for matched filtering, waveforms 
extracted from numerical-relativity simulations are expected to facilitate the detection of 
genuine gravitational waveforms by interferometric detectors. 

Early attempts to evolve the Einstein equations relied on the Arnowitt-Deser-Misner 



(ADM) decomposition |2ll. l22l|. The resulting ADM system proved only weakly hyperbolic 
when expressed in first-order form, a fact partly accounting for difficulties associated with 
its numerical evolution 23j, |24J. Difficulties in evolving black- hole solutions to the Einstein 



equations also stem from singularities, gauge conditions within the computational domain, 
and unstable constraint violation. For over ten years, the goal of accurate and stable numeri- 
cal integration of the Einstein equations has continuously spurred the interest of numericists 
and theorists alike, leading to a wealth of new formalisms j25h52| (this list is not exhaustive). 

To evolve binary black holes, numerical relativists currently use one of the following 




versions of the Einstein equations: the generalized harmonic (GH) system [48j, |49_|, |5l|, [53 



or the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) system 
difference approach with adaptive mesh refinement, Pretorius 



Using a finite- 
used a constraint- 



suppressing second-order form of the GH system (suggested by Gundlach et al. [481 ]) to evolve 



a binary through inspiral, merger and ringdown. Lindblom et al. [5lJ recast the second- 
order GH system into a first-order symmetric-hyperbolic evolution system with constraint 
suppression comparable to that of the second-order system. This first-order GH system 
has been used to successfully simulate binary black holes evolution with nodal spectral 
(pseudospectral) methods [l4 , 15, 54 ] . More recently, Ref. jHEJ has introduced a new penalty 
method for nodal spectral evolutions of spatially second-order wave equations. This work 
provides a foundation for solution of the second-order GH system via spectral methods, 
and has been used to evolve the Kerr solution [56[ and the inspiral of binaries. Typically 
written in a spatially second-order form, the BSSN system [36| has seen widespread use 
by numerical relativity groups that employ finite-difference techniques to evolve binaries. 
Ref. M presented a nodal spectral code to evolve the BSSN system in second-order form. 



58 



The system proved unstable when tested on a single black hole. In more recent work 
longer evolutions were obtained through the adoption of better gauge conditions, filtering 
methods, and more distant outer boundaries. The BSSN system has also been evolved 
in a first-order strongly-hyperbolic formulation for a single black hole with nodal spectral 
methods [H^|. Such evolutions of a single black hole exhibited instabilities similar to those 
reported in Ref. (58| . 

Corresponding to the two versions of the Einstein equations discussed in the last para- 
graph are two distinct techniques for the treatment of singularities in numerical relativity. 
Evolutions based on the GH system have used black-hole excision, whereby the interior of 
an apparent horizon is removed (excised) from the computational domain. This technique 
relies on horizon-tracking and gauge conditions which ensure that inner boundaries of the 
computational domain are pure out-flow, whence no inner boundary conditions are needed. 
Evolutions based on the BSSN system have relied on the moving-punctures technique [!, 



3 



also coined "natural excision." Technically much easier to implement than excision, this tech- 
nique features mild central singularities which evolve freely in the computational domain. 
Initially these puncture points may represent either asymptotically flat regions or "trum- 
pets." Hannam et al. first discussed cylindrical asymptotics in moving puncture evolutions 



60, 61 , see also 6.2-66 



Relative to the alternative systems previously discussed, the BSSN system in second or- 
der form affords an easier treatment of singularities and features a relatively small number 
of geometric variables directly related to the foliation of spacetime into spacelike hypersur- 
faces. However, to date, spectral methods for black-hole binaries have been successfully 
implemented only for the first-order GH system. The binary black hole problem is essen- 
tially a smooth one (singularities reside on sets of measure zero censored by horizons), and 
spectral methods exhibit well-established advantages over finite-difference methods for long- 



time simulation of such problems |67J. Therefore, the development and analysis of a stable 
spectral implementation of the full BSSN system is a worthwhile goal in numerical relativity, 
and the motivation behind the pioneering investigations reported in Refs. [57 59 . 

In Refs. 



50, 62 



Brown introduced a spherically reduced version of the BSSN system as a 
test bed for tractable examination of theoretical and computational issues involved in solving 
this system. Indeed, appealing to the simplicity of this system, he offered geometrical and 
physical insights into the nature of the moving-puncture technique and its finite-difference 
implementation [62, 65, |6~6| (see also in, Here, we exploit this system to a similar end, 
using it as a simplified setting in which to develop spectral methods for the stable integration 
of the BSSN system. Precisely, we develop and test a nodal discontinuous Galerkin method 
(dG) [68| for integration of the spherically reduced BSSN system. While Brown's chief focus 
lay with moving punctures, for further simplicity we adopt the excision technique. Clearly, 
the problem we consider is not as daunting as the one confronted by both Tichy and Mroue 
Nevertheless, our method is robustly stable, and therefore might serve as a stepping 
stone toward a stable dG-based formulation for the full BSSN system. The conclusion offers 
further comments toward this end. 



Nodal dG schemes are both well-suited and well-developed for hyperbolic problems |68 



Although mostly used for hyperbolic problems expressed as first-order systems, dG methods 
have also been applied to systems involving second-order spatial operators, typically via dG 



interior penalty (IP) methods [69H74 



(Refs. [7^-77] discuss the concept of hyperbolicity 
78| | in the context of such systems.) Penalty methods of a different type were exploited in 
Ref. [55[ for the wave equation written in second order form. Local discontinuous Galerkin 
(LDG) schemes, developed initially by Shu and coworkers [79n8l|, constitute an alternate 



approach for integration of spatially second-order systems. LDG schemes feature essentially 
the same auxiliary variables as those appearing in traditional first-order reductions, however 
in LDG schemes such variables are not evolved and arise only as local variables. The basic 
difference between dG IP and LDG methods is the manner in which subdomains are coupled. 
The method we described for the spherically reduced BSSN system is essentially an LDG 
scheme. 

This paper is organized as follows. Section [Til collects the relevant equations from Brown's 
presentation, and develops some further notation useful for expressing the spherically re- 
duced BSSN system in various abstract forms. Section [TTT] presents our nodal dG scheme 
in detail, and Section [IV] documents the results of several numerical simulations testing our 
scheme. Our conclusion discusses possible generalization of our method to the full BSSN 
system. Several appendices collect further technical details. In particular, Appendix [C] 
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considers a simple system which models the spherically reduced BSSN system, giving an 
analytical proof that the model system is L 2 stable in the semi-discrete sense. 



II. SPHERICALLY SYMMETRIC (GENERALIZED) BSSN EQUATIONS 

As shown by Brown [E(j, the BSSN system can be generalized to allow for a conformal 
metric without unit determinant, and this paper focuses on the spherical reduction of this 
system, also considered by Brown in [62|]. In fact, this spherical reduction relies on freedom 
present in the generalized BSSN system, since spherical-polar coordinates should not be 
associated with a unit-determinant conformal metric. Although we work with the spherically 
reduced generalized BSSN system (subject to Brown's Lagrangian condition, to be precise), 
we will nevertheless describe it as the spherically reduced BSSN system. 



A. Basic variables and spherically reduced system. 

The conformal-traceless decomposition of the geometry associated with a spacelike 3- 
surface is 

9ab = X~ l 9ab, K ab = x' 1 (A ab + ^g ab K) , (1) 

where g ab is the physical 3-metric and K ab is the physical extrinsic curvature tensor. The 
BSSN variables are the conformal metric g ab , the conformal factor x, the trace-free extrinsic 
curvature A ab , the trace K = g ab K ab , and the conformal connection T a = —g~ 1 ' 2 d b (g 1 ' 2 g ab ), 
where g is the determinant of the metric. The BSSN system also includes the lapse a, shift 
vector P a , and an auxiliary vector field B a used to define the 'T-driver" for the shift. 

Following Brown, we adopt a spherically symmetric line element, 

ds 2 = -a 2 dt 2 + x~ X 9rr(dr + (5 r dtf + X~ l 9ee(d6 2 + sin 2 6d<j) 2 ), (2) 
along with the spherically symmetric Ansatz: 

/ r r \ / i o o \ 

T a = -cos 9 /{gee sin0) , A ab = A„ -g ee /(2g rr ) . (3a) 

V / \0 -g ee sm 2 6/(2g rr ) J 

Subject to the assumption of spherical symmetry, the basic variables are Xi 9m 9ee, A rr , K, 
r r , a, P r , B r . All are functions of t and r, and satisfy the following spherically symmetric 



5 



(generalized, Lagrangian-form) BSSN systemQ. 

d t a = (3 r a - 2aK - (d t a) (4a) 
d t (3 r = f3 r (3 T ' + - A B r - (d t p r )o (4b) 
d t B r = f3 r B r ' + \{d t T r - f3 r T r ') - rjB r - {d t B r ) (4c) 

d tX = « + \k* x - Pg* - 2 J%* - V x (4d) 

6 6g rr 6gee 6 

d t g r r = g rr + o9rrP - 2A rr a ^ (4e) 

6 6 6g ee 

o 1 or i , A rr g 9S a gggf3 r g' 2 , 

d t goe = ^ g e e + o odooP (4f) 

o g rr og rr o 

RA _ RrA , 4 , W^At 2 j 8V w A T , 2ax(g' rr ) 2 axWee? «(x') 2 
3 3# rr 3^ e0 3^ r 3gj} g 6% 



3 2g rr gee 3g rr 3ggg 6g rr 6ggg 3 3 

2_ xa „ _o^ + oocge_2aAl + ^ _ 2^ (4g) 

3 3g rr 3(?6»6» grr 3(76*0 



^ = /3^ + ^-^ + ^-^ + ^ + ^ (4h) 

2aL g rr gee 2g rr gw ^g„ 3 



Ctl —pi H o h Q 1 3 Q 2 2 

g$ r gee Grrgee g$r 3g rr g l Tr g z rr x 
, 4/? r " W„) 2 | , ggfg ; (4i) 

3prr grr (gee) 2 Q(g r r) 2 3g069rr' 

where the prime stands for partial r-differentiation. Eqs. (HH-i) are Brown's Eqs. (9a-f) 
listed in 1621, subject to his Lagrangian condition (corresponding to v — 1 in Brown's equa- 



tions). The first three equations (HK-c) comprise the gauge sector, and these are essentially 
spherically symmetric versions of the standard "1+log" and 'T-driver" conditions listed in 
Eqs. (1) and (2) of [62|. However, we have introduced the following minor modifications. 
First, (d t a)o designates a constant term which ensures that the right-hand side of the a 
evolution equation (J3^) vanishes at the initial time. This source term as well as the analo- 
gous terms appearing in the evolution equations (JIb,c) for (3 r and B r are needed to enable 
a static evolution of the Schwarzschild solution in Kerr-Schild coordinates. Second, the pa- 
rameter A (perhaps with functional dependence) modifies the hyperbolicity of the first-order 
system j82[. The damping parameter r\ typically appears in standard versions of these gauge 
evolution equations. (See Sections III CI and IIVBI for further discussions.) For this BSSN 
system, we have three constraints: the Hamiltonian constraint "H, the momentum constraint 
Ai r , and the constraint Q r resulting from the definition of the conformal connection T r . In 



1 For this system the determinant g = g r r(gee) 2 sin 4 9 is not unity 
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spherical symmetry, these constraints are written as follows: 



n 

M r 



■2g 



A' 

Qrr 
9rr 



IK 1 3A rr x' 3A rr g0g 



+ 



+ r. 



q 2 

ifrr 



gee g-rr-gee g r 



X9rr9ee 
gl r gee 



x'g'rr , xig, 

q 2 



I \2 



2grrggg 

(5a) 
(5b) 

(5c) 



These expressions are the ones listed by Brown in [62] . Eqs. (Jlfe,f) also ensure that the 
determinant factor g/ sin 4 ^ = g rr {gee) 2 remains fixed throughout an evolution. 



B. Abstract expressions of the system 



We define the following vectors built with system variables: 



u 



( X \ 

grr 

gee 
a 



( B r \ 

A 

K 



( x'\ 



Q 



>rr 



a 



(6) 



Introduction of Q might seem unnecessary at this stage, but proves useful in the construction 
of our discontinuous Galerkin scheme. In terms of the vectors u, v, and Q we further define 



W u:v 



w v:C 



u 



W = W u:Q = I v 

Q 



(7) 



Here we have introduced "colon notation" 83] to represent (sub)vectors and (sub) matrices, 
although we employ the notation over block rather than individual elements. In the first- 
order version of the system (TjJ the components of Q are promoted to independent fields, in 
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which case the corresponding principal part features 
d t B r = f3 r B r " 



3g7 r K + + TuTT^Q^ + w—Qa, 



6(„. 

2 12 

/^rr + ^9rraxT r ' + -OlQ' x ~ ~xQ' a 3g ^g rr 



d t K = (3 r K' 

d t r = p r r r ' 



—Q' 

birr 

AaK' 4 B r B r 



<9 t Q x 

9tQg rr 
dtQaee 



3g rr 3g rr p 6((7 rr ) 2 9rr 
2 ^X Qf 2pXrt 

3g rr 9rr 3g ee 



3gee9rr 



3 



A' - 



3g r 



3gee 

~Qg rr ~~ g 



(8a) 

(8b) 
(8c) 

(8d) 

(8e) 

(8f) 

(8g) 
(8h) 
(8i) 



where all lower-order terms on the right-hand side have been dropped. This sector of prin- 
cipal parts of the first-order system has the form 



d t W v:Q + A(u)W' u:Q = 0, 



(9) 



where (minus) the explicit form of the 9-by-9 matrix A(u) is given below in ( lAlj) . The 
first-order version of (HI) takes the nonconservative form 



d t W + A(u)W = S(W), A(u) 



05x5 


05x9 


09x5 


A(u) 



(10) 



where S(W) is a vector of lower order terms built with all components of W. Partition of 
A(u) = A(u) v: q jV: q into blocks corresponding to the v and Q sectors yields 



(11) 



A{u) vv 






M u )qq 



A{u)- 

Using these blocks, we then define the 9-by-9 matrix 

A(u) = A{u) u , v ^, Q 
and express (j4j) as 

d t W u .. v + A{u)W' v:Q = S{W) 
Q = u', 



05x4 


05x5 


A{u) vv 


M u )vQ 



(12) 



(13a) 
(13b) 



where S(W) = S{W) U . V . 
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field 




V 


/il — u 


-^2,3 


M2.3 = -/3 r 




^4 = ± V 2a X/9rr 


xf 




xi 


^6 = ± V X /drr 



TABLE I. Characteristic speeds. These speeds are the eigenvalues listed in (|A2|) . 



C. Hyperbolicity and characteristic fields 

Although our numerical scheme deals directly with the second-order spatial operators 
appearing in we first consider the hyperbolicity of the corresponding first-order system 
( [TO]) . The characteristic fields and their speeds are found by instantaneously "freezing" 
the fields u in A{u) to some value uq, corresponding to a linearization around a uniform 
state. Below we continue to write u for simplicity with the understanding that u is really 
the background solution Uq. Of primary interest is the range of uq for which the system is 
strongly hyperbolic (75-78 



Appendix [A] shows that the characteristic fields corresponding to (j4]) are as follows: (i) 
all components of u (each with speed 0), and (ii) the fields 

Xx = g ee Q grr + 2g rr Q geg (14a) 

X 2 = g rr T r + \Q X - ^-Q grr - —Q gee (14b) 



ITT 



Xs = 9 ~yB r + \Q X - ^-Q 9rT - ^-Q geg (14c) 
A x l grr gee 

Xf = ±^^K + Q a (14d) 

o j 11 1 

Xt = T—p=Ar ± 1\\—K + 2g rr T r + -Q x Q 9rr + — Q gge (14e) 

yjgrrX \ X X grr gee 



r ± = _^_grr_ B T , a V^gTT K P_ Q 

6 ' 4 A (2a X - A) 8(/^ rr =F vTO ^ 9rr 



with the speeds listed in Table HI To ensure strong hyperbolicity we must necessarily require 

A>0, {/3 r ) 2 g rr -\^0, 2a X -\^0, (15) 
as shown in in Appendix El where further conditions are also given. When A = 1 the 



hyperbolicity condition of Ref. [62| is recovered. In fact, the system could be recast as 
symmetric hyperbolic. Indeed, as it involves one spatial dimension, the relevant symmetrizer 
can be constructed via polar decomposition of the diagonalizing similiarity transformation. 
However, we will not exploit this possibility. 
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This system admits an inner excision boundary provided 

(3 r > max 



1a\ 



<a 2 x 



9r 



9r 




(16) 



holds at the inner boundary. This condition ensures each characteristic field has a nonposi- 
tive speed at the inner boundary, and therefore the inner boundary is an excision boundary 
at which no boundary conditions are needed. The extra flexibility afforded by the param- 
eter A could be used to maintain rigorous hyperbolicity by moving the points at which 
the conditions in f fl5|) are violated outside of the computational domain. Furthermore, for 
A = 1 Eq. fflBT) conceivably fails or is only satisfied close to r = where field gradients are 
prohibitively large. The troublesome X§ gauge mode has a positive speed — j3 r + g rr - 
Indeed, for the conformally flat Kerr-Schild system considered in section IIVBI an inner ex- 
cision boundary is only possible provided A is small enough. 

The transformation (1141) can be inverted in order to express the fundamental fields in 
terms of the characteristic fields: 



B 7 



J 1 , 



1 A 



6 9rr9ee 



r\2 



m 2 9r 



9rrX 

~2a~ 



(*4 



x; 



Xax 



3g rr (2a X - A) 



XI 



(X 5 - x 5 



K 



X 



8ag rr 
1 1 



(X 4 



-X 



4 ) 
r\2 



(17a) 
(17b) 

(17c) 



L(/3 r )Vr-A 

\-xr) 



6 g rr gee 

1 X W r ?9rr ~ 3A 

12 g rr g$e I (/3 r )Vr - A 
X 



1 2 (TV 

Xi + —(x 2 - x 3 ) + - x (x+ + xz) 

g„ 3g rr (2a X -X) 



x 1 + *x 3 -l 7 <** ^ 

2 3(2a X -X) 



+ x 6 

W'f9rr 



Qgee((/3 r ) 2 9rr - A) 



3A 4 2 

X\ + -g rr X 2 — g rr X 3 + 



axg r 



3 (2a X - A) 



(X+ + X 4 



(xt + x;) 



— o9rr(Xs + X 5 ) — -g rr (X£ + X 6 



4#r 



+ 



iP 



r\2 



\2{{^fg rr - A) 



X\ — -geeX 2 



—Qf)nX 3 — — 



aX9ee 
3 (2a X - A) 



(Xt + X7 



Qc 



+ ^g ee (X+ + X 5 ) + -g 6 e(X+ + X 6 
1 



-(X+ + X7) 
(3 r \ 



8grrgee((/3 r ) 2 g rr ~ A) 



A 



(2a X - A) V Sg 



ax 



(X 



x; 




(x e 



x; 



(17d) 
(17e) 
(17f) 

) 

(17g) 
(17h) 

•)■ 

(17i) 
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We will refer to this inverse transformation when discussing outer boundary conditions for 
our numerical simulations in Sec. IIVBI 



III. DISCONTINUOUS GALERKIN METHOD 



This section describes the nodal discontinuous Galerkin method used to numerically solve 
(JIJ). We adopt a method-of-lines strategy, and here describe the relevant semi-discrete scheme 
while leaving the temporal dimension continuous. To approximate (jlj), we follow the general 
procedure first introduced in Ref. 84|. Our approach defines local auxiliary variables Q = u', 
and rewrites the spatially second-order system as the first-order system (fT3"a ). Once we 
use (113b ) to eliminate Q from (fT3a ). we recover the primal equations (j3J). The auxiliary 



variable approach was later generalized and coined the local discontinuous Galkerin (LDG) 
method in Ref. [l9[. We may refer to our particular scheme as an LDG method, but note 
that many variations exist in the literature. We stress that in LDG methods Q is not evolved 
and is introduced primarily to assist in the construction of a stable scheme. 
Equations ( !T2|) and (!T3k) imply that the physical flux function is 



F(W) 



F U (W) \ 
F V {W) ) 



A{u)W v:i 



( 5x i \ 
\f(W)J 



f 



fx 
V fr 



\ 



J 



(18) 



Only the evolution equations for B r , A rr , K, and r r give rise to non-zero components in F, 
and we have collected these non-zero components into a smaller vector / = F v . Inspection 
of OH]) determines these components. For example, from flSfc) we find 

X 



f K = -(3 r K+^Q c 

Qrr 



(19) 



A. Local approximation of the system ([13]) 



Our treatment closely follows 85], but with the equations and notations relevant for 

this paper. Our computational domain Q is the closed r-interval [a, 6]. We cover Q with 

fc miB > 1 non-overlapping intervals D fc = [a h , b h ], where a = a 1 , b = 6 fcmax , and b h ~ l = a k for 
u — o . . . u 

" J i "rnax- 

On each interval D fc , we approximate each component of the system vector W by a local 
interpolating polynomial of degree N. For example, 

N 

X k h (t,r) = J2x(t,r^(r) (20) 

3=0 

approximates x(t,r). Throughout this section, approximations are denoted by a subscript 
h (see j68| for the notation). For example, Wh and fh are approximations of W and /. 
Although Q = u', Qh and u' h are not necessarily the same. In (120]) £j(r) is the jth Lagrange 
polynomial belonging to D h , 

N u 



i+3 
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Evidently, the polynomial Xh interpolates x at r j- To define the nodes rj*, consider the 
mapping from the unit interval [—1,1] to D fc , 

r k (u) = a k + l(l + u)(b k -a k ), (22) 

and the yV+1 Legendre-Gauss-Lobatto (LGL) nodes Uj. The Uj are the roots of the equation 

(l-u 2 )P' N (u) = 0, (23) 

where Pn{u) is the iVth degree Legendre polynomial, and the physical nodes are simply 
r k = r k (uj). In vector notation the approximation f)20p takes the form 

X k h {t,r) = X k {t) T l k {r) 1 (24) 

in terms of the column vectors 

X k (t)= [ X (t } r k ) } --- }X (t,r k N )] T } ^)=K(r),..,4(r)] r (25) 

On each open interval (a h , b k ) C D k and for each component of the equations in ([TBI , we 
define local residuals measuring the extent to which our approximations satisfy the original 
continuum system. Dropping the subdomain label k on the polynomials and focusing on 
the K equation as a representative example, the local residual corresponding to (jUi) is 



\R K )\ = -d t K h + (frK\-C£^ +('- 

V 9rr / h V 



XQg rr Qct\ ( XQggeQa 



h \ 2 9rr J h V 9rr9ee / h 



Q a Qx\ , f3aA 2 \ (I 2 



nmn^;^!' (26) 

Here, for example, the expressions readEI 

{FK>) h = F h Kl =9^9**. (27) 

\ ^9rr J ^9rr,h 

We similarly construct the remaining eight residuals, e.g. (R grr )h and (i?r r )h, as well as five 
residuals corresponding to (|13b ). For example, one of these remaining five is 

(R Qa ) k h = -Q a , h + a' h . (28) 

Let the kth inner product be defined as 

( u > v ) D k = / dru(r)v(r), (29) 

and consider the expression (£ k , (_Rx)^) D fc. We call the requirement that this inner product 
vanish Vj the kth Galerkin condition. For each component of the system and for each k 
there is a corresponding Galerkin condition, in total 9k max (N + 1) equations for (TT3k ) and 



2 At this stage the first expression is generically a polynomial of degree 2N — 1 and the latter is not a 
polynomial. The conventions adopted in Eq. (|27j) prove useful while working with the residual. However, 
later on in Sec. 1111 Cl to obtain the final form (|47[) of the numerical approximation corresponding to (|26[) . 
we will replace nonlinear terms with degree- TV polynomials. 
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5k mSuX (N + 1) for (113b ). Enforcement of the Galerkin conditions on each D k will not recover a 
meaningful global solution, since they provide no mechanism for coupling the local solutions 
on the different intervals. Borrowing from the finite volume toolbox, we achieve coupling 
through integration by parts on r and introduction of the numerical flux f* at the interface 
between sub domains. 

In ( 1261) we only need to consider ((3 r K')h and (xQ' a /9rr)h, as the other terms comprise a 
component of the source vector Sh- Using integration by parts, we write 



(£ k ^K') h ) Qk = - 

{Q,(xQ'j9Tr)h) D >=- 



dr 



dr 



K h \+{F h K h )e 

XhQa,h 



■jk Xh 
9rr,h 



Qa,h 



+ 



9rr,h 



pk 



(30a) 
(30b) 



In these formulas, we have retained the domain index k on £ k , while continuing to suppress it 
on Kh, g rr ,h, etc. Moreover, we have suppressed the r-dependence in all terms on the right- 
hand side. Addition of these formulas along with the definition fx,h = — ((3 r K)h+(xQa/grr)h 
gives 



(£ k j ,(P r K') h -txQ'j9rr) h ) l 



dr 



pk Xh 

■j 

9rr,h 



Qol.I 



— flC.h 



(31) 



In lieu of ( 13T|) . we will instead work with the replacement 
(£ k , ((3 r K') h - ( X Q'j9rr)k) Dk -> - / Or 



{t)fl)'K h - [£) 



')k Xh 



9rr,h 



Qa,h 



* gk 



■ (32) 



This replacement features a component of the numerical flux rather than a component 
fx,h °f the boundary flux. The numerical flux is determined by (as yet not chosen) functional 



/* = f*(W + ,W~ 



(33) 



where, for example, W~ is an interior boundary value [either W k (t,a k ) or Wj^(t,b k )} of the 
approximation defined on D fc , and W + is an exterior boundary value [either W^~ x {t^ b k ~ 1 ) or 
W k+1 (t, a k+1 )] of the approximation defined on either D fc_1 or D k+1 . We discuss our choice 
of numerical flux in the next subsection. We now employ additional integration by parts to 
write the above replacement as 

{£), (f3 r K') h - ( X Q'j9rr)h) Dk -> f dr£ k Uk' - + (f KA - f* K ) £ k (34) 

Ja k \ 9rr J h a 

Rather than the exact fcth Galerkin condition (£ k , {RK)h) Dk — 0, Vj for the K component 
of ( TT3|) on D k , we will instead strive to enforce 



{£ k 3 ARK) k h ) Dk = UK,h-r K )£ k 3 



Vj 



(35) 



3 In the context of the dG method here, + and — denote "exterior" and "interior" , and have no relation to 
the ± using to denote the characteristic fields and speeds in Table |U For characteristic fields and speeds, 
+ and — mean "right-moving" and "left-moving" . 
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although our treatment of nonlinear terms will lead to a slight modification of these equations 
(we return to this issue shortly). The other components of (113b ) are treated similarly, as 
are the components of (fl3b ). Recall that, for example, Q a = a'. Formally using the same 
dG method to solve for Q a , we arrive at the replacement 



dri) 



-Q a ,h + ot h ) - (a h - a* 



(36) 



which again features a component a* of the numerical flux. The auxiliary variables are 
constructed and used at each stage of temporal integration. We then have 



as the corresponding enforced A;th Galerkin condition. 



Vj 



(37) 



B. Numerical Flux 

To further complete our dG scheme we must specify functional forms for the compo- 
nents of the numerical flux introduced in the previous section. We distinguish between the 
physical fluxes (components of /) and the auxiliary fluxes (components of u) arising from 
the definition of the auxiliary variables. These choices are not independent as the resulting 
scheme must be stable and consistent. Our choice follows 
problems. Additional analysis of this flux choice appears in 



86 
68 



which considered diffusion 



87 



Let us first consider the numerical fluxes corresponding to the physical fluxes and of the 
form (|33|) . The numerical flux vector is a function of the system and auxiliary variables 
interior and exterior to a subdomain. A common choice for /* is 

f = {{fh}} + \ W , ^-component of f: f* K = {{f K ,h}} + \ [[K h ]] , (38) 

where, as an example, we have also shown the component of /* corresponding to the analysis 
above. Respectively, the average and jump across the interface are 

{{fh}} = \{f + + D , W = n-tT + n^ + . (39) 

Here r is a position dependent penalty parameter (fixed below) and n~(n + ) is the local 
outward pointing normal to the interior (exterior) subdomain. The role of r is to "penalize" 
(i. e. yield a negative contribution to the L2 energy norm) jumps across an interface. An 
appropriate choice of r will ensure stability, and we now provide some motivation for the 
choice (j4"Tj) of r we make below. 

Were we treating the fully first-order system (fTO]) . the local Lax-Friedrichs flux would 
often be a preferred choice due to its simplicity [68]. In this case, the constant u in the 
numerical flux formula T* = {{J^h}} + I^H^i]] obeys to > max|/i(Viy^ r (W / )) |. Here, 
F{W) = A(u)W, the notation fi(-) indicates the spectral radius of the matrix within, and 
the max is taken over interior W~ and exterior W + states. Motivated by OH]), we adopt a 
similar but simpler prescription, substituting the field gradient 

V Wv:Q A(u)W v:Q = A(u) (40) 



14 



for VwF(W). Precisely, we assume the scaling 

T (b k ) = r{a k+1 ) = T k+1/2 = C ■ max\n(A(u))\, (41) 

where C = 0(1) is a constant chosen for stability. Larger values of C will result in schemes 
with better stability properties, whereas too large a value will impact the CFL condition. 
At the interface point \ k+l l 2 = D fc fl D fc+1 , the vector Uh has two representations: u~ at b k 
and u + at a +1 . The max in (j4"Tl) is taken over the corresponding two sets of field speeds. 
More precisely, the speeds in Tabled are computed for both u~ and u + , and the maximum 
taken over all resulting speeds. For the auxiliary variables, a penalized central flux is used 
The definition with one representative component is 

u* = {{u h }} - - \[u h ]\ , a-component of u*: a* = {{a h }} - - \[a h ]\ , (42) 

with similar expressions for the remaining components. 

We stress the following point. Since the interior coupling between subdomains is achieved 
through the numerical flux forms (1411) and (|42p . the inverse transformation fl 1 Tp expressing 
the fundamental fields in terms of the characteristic fields is not required to achieve this 
coupling. On the other hand, imposition of physical boundary conditions may still rely on 
( I17p . since this transformation allows one to fix only incoming characteristic modes. 



C. Nodal form of the semi-discrete equations 

Let us introduce the kth mass and stiffness matrices, 

M k = f dre k (r)£ k (r), 5* = f drt k {r)^{r). (43) 

These matrices belong to D fc , and the corresponding matrices defined on the reference interval 
[—1,1] are 

Ma = J duti(u)tj(u), Sij = J duti(u)^(u), (44) 

where £j(u) is the jth Lagrange polynomial determined by the LGL nodes Uj on [—1,1]. 
These matrices are related by M k - = \{b k — a k )Mij and S k j = S^, whence only the reference 
matrices require computation and storage. 

We will use the matrices M k and S k in obtaining an ODE system from ( J26l) and (!3~5T) . 
Towards this end, we first approximate the nonlinear terms (products and quotients) in 
(l2"u]) by degree-iV interpolating polynomials. Such approximations are achieved through 
pointwise representations. For example, (Q a Q x / g rr )h appears in (1261) . and is expressed in 
the following way: [cf. footnote |2] 

(QaQx\ t. \ (Qa,hQ X ,h\(. \ v Qa,h if, T k ) Q %h (t, T k ) , 

VT^) ( ' r)= [~tt~) { ' r) Q h (tr k ) ^ (r) ' (45) 



Note that the expressions on the right and left are not equivalent due to aliasing error |67 
Our vector notation for this replacement will be 

QaQx\ (j. __\ . ( QaQx \ ,.,././' (>!< 
\ 9r 



(t,r)^ P^ ){tyt{r). (46) 
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Operations among bold variables are always performed pointwise. Making similar replace- 
ments for all terms in (1261) . and then carrying out the integrations in ( 135]) . which bring in 
M k and S k , we arrive at 



d t K = (3 r DK 



X DQ a , 1 XQg rr Qa XQ gee Qc 



+ - 



lQ a Q 



Qrr 

3 aAl. 



q 2 

Urr 



9 rr 9ee 



2 9r 



2 9 



+ l - (X K 2 + M' i i k {f K , h -r K ) bk 



(47) 



where we have again suppressed the superscript k on all terms except £ h (r), and the sub- 
script h is dropped on all boldfaced variables. As described in 68|, the spectral collocation 
derivative matrix 



dr 



(48) 



can also be expressed as D k = (M k )~ 1 S k , which appears above. Eight other semi-discrete 
evolution equations are similarly obtained, with nine in total (one for each component of 
W u:v ). Additionally, we have 



Da + M~ 1 £ k (a* - a 



(49) 



as one of the auxiliary equations, with five in total (one for each component of Q = Wq-q) 



D. Filtering 

Like other nodal (pseudospectral) methods, our scheme may suffer from instabilities 
driven by aliasing error j67[. Filtering is a simple yet robust remedy. To filter a solu- 
tion component, such as x> we use the modal (as opposed to nodal) representation of the 
solution: 

N N 

X k h (t,r) = Y,x(t,r k )£ k (r) = £xj(*)^to, (50) 

3=0 3=0 

where Pj(r) is the jth Legendre polynomial. Let rjj = j/N, and define the filter function 

[1 for < f]j < NJN 

^ H \exp(-e(^|f) 2 ') for NJN < r» < 1. (51) 

At each timestep we modify our solution component according to 

N 

A ^ W )" = ^ a(%) ^ (i) p. (r) . (52) 

3=0 

Evidently, the modification only affects the top N — N c modes, and is sufficient to control 
the type of weak instability driven by aliasing [68j . The numerical parameters N c and e are 
problem dependent. For our simulations, we have taken e ~ — log(e mac h) = 36, where £ ma ch 
is machine accuracy in double precision. 
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E. Model system 

To better illustrate the basic properties of our method, we consider a toy model. Namely, 
the following spatially second-order system: 

d t u = u + av — u 3 + g(t, x) (53a) 
d t v = u" + v' - (u + v)(u) 2 + v 2 u 2 + h{t,x), (53b) 

where a > 1 is constant and g and h are analytic source terms to be specified. In contrast 
to (J6]), here u, v, and Q = u' are scalars rather than vectors. System (1531) admits a first- 
order reduction in which u' is defined as an extra variable. Since this first-order reduction 
is strongly hyperbolic, the spatially second-order system (1531) is also strongly hyperbolic by 
one of the definitions considered in [77]. The characteristic fields X ± and speeds fi^ are 

X + = yfav — u, /i + = \fa — 1; X~ = \fav + u', ji~ = —{\fa + 1). (54) 

To construct a local dG scheme for this system, we first rewrite it as 

d t u = Q + av — u 3 + g(t, x) (55a) 
d tV = Q' + v' - {u + v)Q 2 + v 2 u 2 + h{t, x) (55b) 
Q = u . (55c) 

Evidently, / = — (Q + v) is the v -component of the physical flux vector 

nv,Q)=(Z U ) = (°V (56) 



FvJ \f 

Note that F has the same structure as (u, v) T . Borrowing from the presentation for the 
BSSN system, we write the analogous semidiscrete scheme on each subdomain D k for the 
model system: 

d t u = Q + av - u 3 + g(t) (57a) 

fok 

d t v = DQ + Dv-(u + v)Q 2 + v 2 u 2 + h(t) + M- 1 £ k (f h -f*) , (57b) 



Q = Du + M- l t{u* -u h ) 



(57c) 



Here, we have suppressed the subinterval label k from all variables except for the vector 
£ k of Lagrange polynomial values. Moreover, following the guidelines discussed above, the 
numerical fluxes are given by 

r = {{h}} + N] - u * = {{«*» - \ [KI • (ss) 

Appendix O analyzes the stability of our scheme, for a more general numerical flux choice, 
as applied to ( 1531) with the nonlinear and source terms dropped. 



IV. RESULTS FROM NUMERICAL SIMULATIONS 



This section presents results found by numerically solving both the model system ( 153]) 
and BSSN system (j3J) with the dG scheme presented in Sec. IIHI 
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FIG. 1. Spectral convergence of fields for model PDE. Respectively, for N = 3,6,9,12, 
a timestep of At = 0.0578, 0.0178, 0.0084, 0.0049 has been chosen for stability and accuracy. In the 
title headings, for example, Ait = u numer — tt eX act- 



A. Simulations of the model system 

The semi-discrete scheme (157"]) has been integrated with the classical fourth-order Runge- 
Kutta method. When integrating this system, we have first constructed Q at each Runge- 
Kutta stage, and then substituted into the evolution equations (joTh.b) for u and v. The 
problem has been solved on a computational domain [0, 47r] comprised of two subdomains 
with a timestep chosen small enough for stability. The initial data has been taken from the 
following exact solution to (153"]) : 



,(t,x) 

^exact {tj 

g(t,x) 
h(t, x) 



[sin (a; — ji t) — sin(a; — fx + t 



[ sin(x — \i t) + sin(x — fi + t 



2v^ 



u. 



exact 



("^exact + ^exactX^exact) V exact U 



exact ) 



(59a) 
(59b) 

(59c) 
(59d) 



where the speeds jj, are found in (1541 . Specification of the boundary condition at a physical 
endpoint amounts to choosing the external state for at the endpoint. We have considered 
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FIG. 2. Scaling of maximum stable At with N for model PDE. 



two possibilities: (i) the analytic state (Q + ,v + ) = (Q exact, v exac t) and (ii) an upwind state. 
For example, at x = An the upwind state igj 

Q Qupwind [("^ )cxact (-^ )numer] ; ^ ^upwind ~ ^ [C^ )exact "I - (-^ )numer] • 

(60) 

Either choice of (Q + , u + ) leads to similar results, and the plots here correspond to the 
analytic state. Figure [T] clearly shows spectral convergence with increasing polynomial order 
N across all fields for the case a = 2. Other values of a, including a = 1 for which X + 
is a static characteristic field, have also been considered with similar results. Appendix O 
demonstrates that our proposed scheme for the system ( 1571) with nonlinear and source terms 
dropped is stable in a semi-discrete sense. Nevertheless, the fully discrete scheme, obtained 
via temporal discretization by the fourth-order Runge-Kutta method, is still subject to the 
standard absolute stability requirement. Namely, if is any eigenvalue corresponding to 
the (linearized) discrete spatial operator, then a necessary condition for stability is that 
/i/jAt lies in absolute stability region for Runge-Kutta 4. We here show empirically that the 
associated timestep restriction scales like N~ 2 , i.e. At = 0(N~ 2 ) for stability. We note that 
such scaling is welcome in light of the second-order spatial operators which appear in the 
system, and suggest a possible worse scaling like A r ~ 4 . Fig. [2] plots the maximum stable 
timestep for a range of N, demonstrating the iV -2 scaling, in line with behavior known from 
analysis of first-order systems (68[. This scaling also holds for the BSSN system. 



4 We remind the reader that, unfortunately, the ± on means something different than the ± indicating 
exterior/interior dG states [cf. footnote |3]. 
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B. Simulations of the BSSN system 

This subsection documents results for simulations of the unit-mass-parameter (M = 1) 
Schwarzschild solution ( 1B9I) expressed in terms of ingoing Kerr-Schild coordinates. Since the 
solution is stationary, temporal integration of the semi-discrete scheme has been carried out 
with the forward Euler method which the dissipation in our method allows. The r-coordinate 
domain [0.4, 3.4] has been split into 3 equally spaced subdomains, and we have set rj = 10, 
A = 0.1, and C = 2 [cf. Eq. (I4"TT) ]. For all simulations At has been chosen for stability. With 
the chosen A, the inner physical boundary r min = 0.4 is an excision surface. At each timestep 
we have applied an (order 2s = 20) exponential filter on the top two-thirds of the modal 
coefficient set for all fields except for g rr and gg$. For stability, we have empirically observed 
that g rr and gee must not be filtered. A detailed understanding of this is still lacking. 

Issues related to physical boundary conditions are similar to the one encountered in 
Sec. IIV Al for the model problem. Similar to before, we have retained Eqs. (I38f42p as the 
choice of numerical flux even at the endpoints. Therefore, at an endpoint the specification 
of the boundary condition amounts to the choice W + of external state. We have typically 
chosen the inner boundary of the radial domain as an excision boundary, and in this case 
W + = W~ is enforced at the inner physical boundary. At the outer physical boundary, 
for W + we have again considered two choices: (i) W^exact and (ii) W^upwind- To enforce 
choice (ii) the inverse transformation (Tl7|) must be used with incoming characteristic fields 
fixed to their exact values, similar to f )60|) . We have tried various versions of choice (ii), 
and in all cases the resulting simulations have been unstable. We therefore present results 
corresponding to choice (i). Although the choice of an analytical external state W ex&ct 
at the outer boundary is stable for our problem, such a boundary condition is unlikely 
to generalize to more complicated scenarios involving dynamical fields. Indeed, the issue 
of outer boundary conditions for the BSSN system is an active area of research, with a 
proper treatment requiring fixation of incoming radiation, control of the constraints, and 



specification of gauge (see Ref. |88[ for a recent analysis). 

For BSSN simulations, our main diagnostic is to monitor the Hamiltonian, momentum, 
and conformal connection constraints. Figure [3] depicts long-time histories of constraint 
violations, whereas Figs. H] and [5] depict long-time error histories for the individual BSSN 
field components. From the middle plot in Fig. El we infer that, up to the indicated numerical 
error, the factor gj sin 4 # = g r r(gee) 2 remains at its initial fixed profile r 4 throughout the 
evolution. These figures indicate that the proposed scheme is stable for long times, and 
exhibits spectral converge with increased polynomial order N. Similar results are recovered 
from M = Minkowski initial data. The stability documented in these plots does not appear 
to rely on inordinate parameter tuning. For example, with the fixed parameters described 
above, we obtain similar plots if we individually vary (i) r min over {0.325, 0.35, 0.4, 0.475} 
(values still corresponding to an excision surface for the given choice of A), (ii) n over 
{1, 3, 7, 10}, (iii) s over {8, 9, 10}. With the polynomial order N ranging over {23, 26, 29, 31}, 
both stability and qualitatively similar exponential convergence is achieved with a single 
subdomain. Likewise, adoption of a larger coordinate domain with more subdomains does 
not significantly impact our results. However, for much larger r max stability requires a 
smaller time step or a time stepper better suited for wave problems (e.g. Runge Kutta 4). 
Finally, we have considered the addition of random noise to all field components at the 
initial time. Precisely, with the system component x as an example, we have set 



(61) 
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FIG. 3. Spectral convergence of constraint violations for M = 1 Kerr-Schild 
initial data. Respectively, for N = 11, 14, 17, 19, a timestep of At ~ 0.0041, 0.0026, 0.0018, 0.0013 
has been chosen for stability and accuracy. 

where each component (nodal value) of Sx{0) is 10 -5 times a random variable drawn from a 
standard normal distribution. Such perturbed initial data also gives rise to stable evolutions. 



V. CONCLUSION 

We have introduced a discontinuous Galerkin method for solving the spherically reduced 
BSSN system with second-order spatial operators. Our scheme shares similarities with 
other discontinuous Galerkin methods that use local auxiliary variables to handle high-order 



spatial derivatives (68|, 79-81, 84, [8^, 87||, and which have typically been applied to either 
elliptic, parabolic, or mixed type problems. The key ingredient of a stable dG scheme is an 
appropriate choice of numerical flux, and our particular choice has been motivated by the 
analysis presented in Appendix O When used to evolve the Schwarzschild solution in Kerr- 
Schild coordinates, our numerical implementation of the BSSN system PJ is robustly stable 
and converges to the analytic solution exponentially with increased polynomial order. By 
approximating the spatially second-order form of the BSSN system, we have not introduced 
extra fields which are evolved. Evolved auxiliary fields result in new constraints which may 
spoil stability. Our main goal has been stable evolution of the spherically reduced BSSN 
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system as a first step towards understanding how a discontinuous Galerkin method might be 
applied to the full BSSN system. Towards that goal, we now discuss treatment of singularities 
and generalization of the described dG method to higher space dimension. 

To deal with the fixed Schwarzschild singularity, we have used excision which is easy in 
the context of the spherically reduced BSSN system. However, excision for the binary black 
hole problem in full general relativity requires attention to the technical challenge of hori- 
zon tracking. State-of-the-art BSSN codes avoid such complication, relying instead on the 
moving-puncture technique. While the moving-puncture technique does involve mild central 
singularities, it may still prove amenable to spectral methods. Indeed, spectral methods for 
non-smooth problems is well-developed in both theory and for complex applications. Since 
the moving-puncture technique can be performed in spherical symmetry [62J, a first-step 
toward a spectral moving-puncture code would be to implement a moving puncture with 
the nodal dG method described here. Such an implementation may adopt Legendre-Gauss- 
Radau nodes on the innermost subdomain, thereby ensuring that the physical singularity 
does not lie on a nodal point (in much the same way finite difference codes use a staggered 
grid). Beyond traditional excision and moving punctures, one might construct smooth ini- 
tial data via the turducken approach to singularities. However, in combination with 1+log 
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FIG. 5. Spectral convergence of solution violations for M = 1 Kerr-Schild initial 
data. See the caption of Fig. U]for details. 



slicing and the Gamma-driver shift condition, turduckened initial data will evolve towards 
a "trumpet" geometry [ 



2jJ. 

Discontinuous Galerkin methods for hyperbolic problems in two and three space dimen- 
sions are well-developed. A generalization of the method described here to three-dimensions 
and the full BSSN system would likely rely on an unstructured mesh. Appropriate local poly- 
nomial expansions for the subelements are well-understood, as are choices for the numerical 
fluxes which would now live on two-dimensional faces rather than single points. Whether or 
not it would ultimately prove successful, generalization of our dG method to a higher dimen- 
sion would rely on an established conceptual framework. Further computational advances 
of relevance to a generalization of our dG method to the full BSSN system (possibly includ- 
ing matter) may include mesh /ip-adaptivity, local timestepping, shock capturing and slope 
limiting techniques [68]. Moreover, recent work [9l[ indicates that enhanced performance 
would be expected were our scheme implemented on graphics processor units. 
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Appendix A: Hyperbolicity of the first-order system. 



This appendix analyzes the matrix A(u) appearing in ( 1101) in order to construct the 
characteristic fields f[T4|) . In matrix form the sector of the principal part of (fTUj) reads as 
follows: 
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(Al) 

which defines the matrix A(u) appearing in (Q, and so also the matrix A(u) in ({TO]) . Note 
that in the last equation the matrix within the square brackets is — A(u). For certain 
configurations of u and A, the system ( fiTJl) is strongly hyperbolic [78j |. that is A{u) has 
a complete set of eigenvectors and real eigenvalues. Indeed, five eigenpairs of A(u) are 
trivially recovered upon inspection of A(u)'s leading 5x5 diagonal block. These correspond 
to eigenvalue and the left eigenspace {£j = ej : 1 < j < 5}, where tj are the canonical 
basis vectors. Since each component of u arises as eJW, each is also a characteristic field. 

The remaining nine eigenpairs are determined by A{u). The eigenvalues of A(u) are 

Ml = 0, /i2,3 



/i4 



Qrr V Qrr V Qrr 



(A2) 



24 



and the corresponding left eigenvectors are 



Xl = (0,0,0,0,0,^,2^,0,0) (A3a) 
x 2 = ( 0, 0, 0, g rr , -, -— , 0, ) (A3b) 

x 3 = 0, 0, 0, -, -— , 0, o) (A3c) 



x£ = (o,0,±y^^,0,0,0,0,l,oj (A3d) 
4 = ( 0, =F— J=, ±2,/^, 2<? rr , -, -— , — , 0, o) (A3e) 
^=14^0^^0,0, £ 



4 A ' ' (2« X -A)' ' ' SCS^TaA^)' 

(!V ,±*/^Y (A3f) 



4gee(P r grr =F v 7 ^) ' ( 2 «X - A) ' V A 

where for example x^A(m) = fM 5 x 5 . Assuming that g rr , ge9, Xi an d a are everywhere strictly 
positive, the eigenvalues are real and the eigenvectors are linearly independent provided 
that (IT51) holds. These eigenvectors are easily extended to eigenvectors of A(u), e. g. as 
Xq — > (0i X 5,Xg ). Then, for example, the characteristic field 



X± = (0 lx5 ,4)W = xtW vlQ , (A4) 



and similarly X? = XjW v: q for j — 4,5 and X k = XkW v: Q for k — 1, 2, 3. The characteristic 
speeds for these fields are and With this convention the speeds listed in Table fl] 
correspond to the Xk and X- in (ITI]) . 



Appendix B: Schwarzschild solution in conformal Kerr-Schild coordinates. 

In Kerr-Schild coordinates, here the system directly related to incoming Eddington- 
Finkelstein null coordinates, the line element for the Schwarzschild solution reads 

ds 2 = -a 2 dt 2 + (1 + 2M/R)(dR + /3 R dt) 2 + R 2 d8 2 + R 2 sin 2 6d(f) 2 , (Bl) 

where R is the area radius, a = (1 + 2M/R)- 1 ' 2 is the lapse, and R = 2M/(R + 2M) is 
the shift vector. The physical spatial metric g a b is the spatial part of this line element. 

To define the corresponding solution to the BSSN system, we use equation g ab = X9ab to 
define the following relationship between line elements: 

dr 2 + r 2 {d6 2 + sin 2 6d<f) 2 ) = + 2M/R)dR 2 + R 2 dQ 2 + R 2 sin 2 6d<p\ (B2) 

so that 

4 + f)(f) 2 ^ 

Then we have 

/ 2M\ l/2 dR dr . . 
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with integration yielding 



r = ^il + ] /l + ^\ e ^V^m } (B5) 

where the constant of integration has been chosen so that the R, r — > oo limits are consistent. 
The second relation in (1B3|) shows that 



2 



1 / / 2M\ 4 

X 



.. 1 + — e 4-4yiT2M7i? -4 = (B6 ) 

16 ^ V R ) X l + y/l + 2M/R K } 

The extrinsic curvature tensor is specified by the expression for K given in flB9h ). the identity 
K = K% + 2K e e , and 

2MV 1/2 2M 



2M\~ 1/2 / R + M \ 2M 



Since = we compute that 

A T = A' - 2K» = + ) (B8) 

Next, since K rr = g Tr K r r = x~ l K^ we have K r r = A rr + ^g rr K. This implies A rr = K'^. — ^K, 
from which we get 0B9g ). In all we have 



2M 



-1/2 



a = (1 + — \ (B9a) 

9rr = 1 (B9C) 

g ee = r 2 = X R 2 (B9d) 



1 / 2M 



X 



4 



16 + V 1 + IT J e4 " 4v ^^ ( B9e ) 

fi r = (B9f ) 

/ 2M\" 1/2 AM f2R + 3M\ , n . 

2MV 3/2 / 3M\ 2M 



A ' = l 1 + irJ 1 1 + itJip (B9h) 

!- = -; = -^s- («») 



To differentiate these expressions with respect to r, we use the identity 

d_R /2 / 2A 

rfr X \ R 

along with the chain rule. 



dR 2M\ _1/2 

~iz = x~ 1/2 l + -5- (Bio) 
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Appendix C: Stability of the model system 



The following stability analysis for the model system f )53|) has been inspired by |80l . 181 . 
After dropping all nonlinear source terms, the system ( |53l) becomes 

d t u = u' + av (Cla) 
d t v = u" + v'. (Clb) 

This section analyzes the stability of (IClj) . considering both the continuum system itself as 
well as its semi-discrete dG approximation. The latter analysis offers some insight into the 
empirically observed stability of our dG scheme for the spherically reduced BSSN equations. 



1. Analysis for a single interval 

Throughout we work with the L2 inner product and norm, 



{f,g)o= / fg, II/IId = VifJh, (C2) 



where D is the spatial coordinate interval (here D may represent a subdomain D k or the 
whole domain Q), and we have suppressed all integration measures. For the continuum 
model we will establish the following estimate: 

\\u'(T, -)|| D + a\\v(T, < C(T) (\\u'(0, -)|| 2 D + a\\v(0, , (C3) 

where the time- dependent constant C(T) is determined solely by the choice of boundary 
conditions. To show (IC3j) . we first change variables with v = i/af, thereby rewriting (IClj) 
in the following symmetric form: 

d t u = u + y/av (C4a) 

d t v = vW' + v'. (C4b) 

Equations ( 1C4k .b) then imply 

l -d t [ {u' f= [ u'{u" + ^v)= [ ^u'v' + ~ [ {u' f (C5a) 

2 JD JD JD * JdD 



2 Jo JD JD * JdD 



(C5b) 



Here vv' and u'u" have been expressed as exact derivatives and then integrated to boundary 
terms, the second equation employs an extra integration by parts, and with only one space 
dimension J dD denotes a difference of endpoint evaluations. Addition of Eqs. ( 1C5k .b) gives 

U [ [v 2 + (u>) 2 } =U [v 2 + (u) 2 + 2^au'v] . (C6) 

^ JD 1 JdD 

Substitutions with the identities 

[{) 2 + [u') 2 ] = \[{v + u) 2 + (v- u') 2 ] , 2u'v = l[(v + u') 2 -(v- u') 2 } (C7) 
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and replacements to recover the original variable v = vj \fa yield 



9D 



at; - «') • (C8) 



9D 



From (IC8p we deduce that the time-dependent constant C(T) in (1C3P must satisfy 



1 + 



av — u 



A2 



|u'(0,-)|| 2 + a|K0, 



< C(T). 



(C9) 



For periodic boundary conditions, we may choose C(T) = 1. Moreover, if a > 1 and 
«' = — y/av is specified at <9D + , then \\u'(t, -)|| 2 + a||f (t, -)|| 2 decays. 

Still working on a single interval (subdomain), we now consider the semi-discrete scheme 
for ( 1C4I) . i. e. ( 1571) with all nonlinear source terms dropped, and with v replaced by vjyfa. 
Derivation of a formula analogous to dQ8|) is our first step toward establishing L 2 stability 
of the semi-discrete scheme. While (|5T|) features vectors, for example u(t), taking values at 
the Legendre-Gauss-Lobatto nodal points, here we work with the numerical solution as a 
polynomial, for example Uh{t,x). These two representations are related by the Lagrange in- 
terpolating polynomials for the nodal set, here taken to span both the space of test functions 
and the space of basis functions. Our scheme is 



Tpd t u h = / ip(Q h + \fav h ) 



£d t v h 



JD k JdD k 



¥>Qh 



<pu'h + (p(u*- U h ) 



(ClOa) 
(ClOb) 
(ClOc) 



where ip, and (p are polynomial test functions. These test functions are arbitrary, except 
that they must be degree- N polynomials. In ( 1C10I) the variables Uh, Vh an d Qh should also 
carry a superscript k, but we have suppressed this. Derivation of a formula analogous to 
( 1C8j) is complicated by the fact that Qh is not evolved. Nevertheless, at a given instant t we 
can assemble Qh from (IClOb ). 

Mimicking the calculation flC5b ) from the continuum case, we first use ( IClOb ) with £ = 
to write 



£ Ink 



'D k & JdD 

The right-hand side of flC5a ) is analogous to 



(VaQ h + v h )v' h + / (y/EQ* + v*)v h 

D JdD k 

V^QhV h + \ f [2(V^Q* + v*)v h - v 2 h ] . 

£ .lank 



(Cll) 



k/ Ql 

D k 



QhdtQh- 



D k 



(C12) 



However, since Qh is not evolved, the term d t Qh must be given a suitable interpretation. 
On the right hand side of (IClOb ) only Uh, u' h , and u* necessarily depend on time, since the 
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test function ip need not be time-dependent. Furthermore, u* is explicitly given as a linear 
combination of Uh, as seen in Eq. (1C21cj) below. Choosing ip = £j, taking the time derivative 
of (IClOb ). and appealing to the commutivity of mixed partial derivatives, we therefore arrive 
at 

/ iAQh= I t j (d t u h y+ [ ^((dtuy-dtUh), (ci3) 

JD k JD k JdD k 

where (d t u)* depends on d t Uh in precisely the same way that u* depends on Uh- We have 
written £j rather than ip in the last equation to emphasize that the result also holds for any 
linear combination of £j (for example <p), and even for time- dependent combinations. Since 
Qh is itself such a combination, we obtain 

\dt! Ql= f Qh(d t u h )'+ f ((d t uy - d t u h )Q h 

^ JD k JD k JdD k 

= [ Qh(Q' h + V^v' h )+ [ {(d t u)* - d t u h )Q h (C14) 

JD k JdD k 

= [ V^QhV h + l [ [2((d tU y-d t u h )Q h + Q 2 h ]. 

JD k Z JdD k 

Addition of flCTTj) and flCTij) gives 



\d t [ (Ql + vl) = \ I [Ql -v 2 h + 2(v^Q* + v*)v h + 2((d tU y - d t u h )Q h ] , (C15) 



the aforementioned analog of (IC8I) . This formula holds on a single subdomain D , and we 
now combine multiple copies of it, one for each value of k. 



2. Analysis for multiple intervals 

To facilitate combination of ( 1C15I) over all k, we change notation. At every subdomain 
interface \ k+l / 2 = dD k R dD k+1 , let the superscripts L and R denote field values respectively 
taken from the left and right. Then the fields evaluated at \ k+l / 2 which belong to D k will be 
M fc+i/ 2 ; and Qk+1/2, while those belonging to D k+1 will be u* +1/2 , v£ +1/2 , and Q§ +1/2 . 

However, at | fc_1 / 2 the values taken from D fe are u^_ ± ^ 2 , v£_y 2 , and Qf_ lj/2 - Note that we 
have also replaced the subscript h, denoting a numerical solution, with k ± 1/2, denoting 
the location of the endpoint value of the numerical solution. With this notation, we define 

= \ [{Q L a f - [v L a f] + (v^Q; + vl) v L a + [(d t u a y - d t u L a ] Q L a , (C16) 

and similarly for A^. The same numerical fluxes appear in both and A^ (i.e. each 
numerical flux takes the same value on either side of an interface), whence fluxes like Q* a do 
not carry an L or R superscript. In terms of these definitions (1015|) becomes 

\dt I (Ql + v 2 h ) = A L k+1/2 -Al 1/2 . (C17) 

L JD k 
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Summation over all D k yields 
2 



-i "•max p "-max 

oft E / + **) = E ( A W> - a £h/J + A L x+ i /2 - Af /2 

fc=l >/U k=l 

^max 1 

= E ( A f- A fc)Ui/ a + A L, + i/2-Af /2 . (C18) 



fe=l 



We have reverted to /i-notation denoting the numerical solution, since the L, R superscripts 
indicate unambiguously the relevant domain used for evaluation at \ k+l l 2 . 
We again seek an estimate of the form 

^max ^max 

X) OIIS* + a||t; fc (r, -)IIdO < ^(=0 X) (H^C * OIId* + all«*(0, OIId*), (C19) 

k=l k=l 

that is essentially the same as the one f)03[) considered in the continuum case. Assume that 
the chosen boundary conditions ensure A^^ +1 ^ 2 — Aw 2 is bounded by a time-dependent 
constant which does not depend on the numerical parameters N and h (subdomain width). 
Establishment of stability then amounts to showing that the remaining sum over interface 
terms in (10181) is non-positive; whence this remaining sum is consistent with C(T) < 1, 
although the boundary conditions may give rise to a different bound. In fact, we will choose 
the numerical fluxes such that each individual interface term is non-positive. At interface 
|fc+i/2 an( j j n ^ ft no t a tion, the jump and average of Vh, for example, are 

I (£+ + d -) = {{dh}} = I (^ +i/2 + d R +i/2 ) ( C 20a) 
n-v- + n + v+ = [[v h ]] = v L k+1/2 - v* +1/2 . (C20b) 



Consider numerical fluxes of the form 



Q* = {{M-|M (C21a) 



v* = {{v h }}-j[M (C21b) 
{M-yW (C21c) 

(d t u)* = {{d t U h }} - y ptU h ]} , (C21d) 

where (jC21b ) induces ( 1C21d ) and where the penalty parameters r u , r v , and tq are real 
numbers. The fluxes defined in fl58l) correspond to r u — 1, r v — 1 + y/a, and tq = 0. In 
terms of these quantities the kth interface contribution in (IC18I) is 

- asi,^ = \ ([[Qi]] - kid + amm - T iiM 2 

1 1 (C22) 

+ y/E{{Q h }} [[vh]} ~ ^ M [M] - {{Qh}} [[d t u h ]\ - | M [[Q h ]} , 

where we have suppressed the k dependence of the right-hand side. Now consider the term 
[[<9tW/i]] • Because d t Uh and Qh + \/adh are both polynomials of degree N, Eq. ( IClOa ) implies 



u 
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FIG. 6. Stable evolutions for the model system. For fixed r v = 1CT 6 and t v = 1 + y/2 
respectively, the left and right plots depict stable choices (determined empirically) of t u and tq 
for the linear model system (|Clh . The stable regions are colored black, but the jagged edges result 
from the discretization of the (t u , TQ)-plane. 



the vector equation dtu = Q + y/av, that is pointwise equivalence on the nodal points of 
D k , which in turn implies [[^tw^]] = [[Qh + v 7 *^]] • Upon substituting this identity into the 
last equation, we arrive at an expression which features only dh and Qh, 

K - A*)| |fc+1/2 = \ ([[Ql]} - [K]]) + {{v h }}[M - T iiM 2 

+ MiQh}} [M - M [[v h ]] - {{Qh}} [[Qh + V^v h ]] - y [[Qh + v^>ft] [[Qh]] ■ 

(C23) 

The i dentities {{^}}[[^]] = \ [[vl]\ and [[Q h + y/av h ]] = [[Qh]] + \/a[[vh]] then simplify 
f lC23|) to 



(AJ - Af)| |fc+1/2 = -| [[v h ]] 2 - MTu 2 +TQ) [[Qh]] [[h]] ~ f [M] 2 . (C24) 

The role of a penalty parameter is now clear. Positive values of r v penalize jumps in dh 
through a negative contribution to the energy. Likewise, positive values of r u penalize 
jumps in Qh through a negative contribution to the energy. However, because the sign of 
[[Qh\] [[vh\] can be positive or negative, only the choice tq = — t u yields an expression for 
(A^ — A^)|[fc+i/2 which is manifestly negative for t u > and r v > 0. A simple estimate 
based on Young's inequality with e (that is, 2a(3 < e^ x a 2 + e/3 2 , where a, (3 > and e > 0) 
shows that for tq = the choice r v > ar u /A also yields stability. 

Figure [6] depicts certain choices of stable penalty parameters for the linear model system 
evolved to tfi na i = 1000 (with a = 2, iV = 10, and At ~ 0.0553), as determined empirically 
with simulations similar to those described in Sec. IIV Al The left plot corresponds to a 
small r v = 10~ 6 , for which the choice r u — 1, tq — is not stable, as expected from the 
theoretical analysis. However, the right plot corresponds to r v — 1 + \fa, for which r u — 1, 
tq = is stable. Motivated by the numerical flux choices f l38f42]) used for the BSSN system 



31 



01]), we have (as mentioned above) set r u = 1, r v = 1 + a/o, and rg = in simulations of 
the nonlinear model ( |53|) . For the nonlinear model system ( 153|) . the theoretically motivated 
choice tq = —t u also yields numerically stable evolutions when t u > and r v > 0. 

For the nonlinear systems (jlj) and (1531) . we do not attempt a formal stability proof. 
Nevertheless, the results of this appendix have served as a guide for our choices of penalty 
parameters. For the BSSN system (jlj), u, v, and Q are block indices [cf. Eq. (JBJ)]. Similar 
to the model problem, we have penalized Q with t u = 1, with chosen large enough to 
heuristically overcome the cross-terms of indefinite size that arise from tq = (we interpret 
equations like t u = 1 componentwise). An analogous choice "tq = — r u " for the BSSN 
system might be possible, but would be considerably more complicated. Indeed, such a 
choice likely entails a matrix of penalty parameters, but we do not give the details here. 
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